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Abstract 

A  new  cohesive  zone  model  is  developed  to  simulate  the  ductile  to  brittle  failure  transition  in  polycarbonate.  The 
model  is  formulated  so  that  as  rate  or  stress  state  changes  within  a  simulation,  the  fracture  energy  and  thus  fracture 
mode  may  also  change  appropriately.  The  ductile  to  brittle  transition  occurs  when  the  cohesive  opening  rate  is  greater 
than  a  threshold  opening  rate  and  when  the  stress  state  is  close  to  plane  strain  in  a  fracture  specimen.  These  effects 
are  coupled  in  a  phenomenological  model.  The  principal  contribution  of  this  work  is  that,  in  a  glassy  polymer,  the 
transition  from  slow  to  fast  crack  growth  as  test  loading  rate  and/or  sample  thickness  are  varied  can  be  simulated 
by  using  a  single  set  of  bulk  and  cohesive  zone  material  parameters.  This  result  enlists  the  use  of  the  Simplified 
Potential  Energy  Clock  (SPEC)  model  and  the  new  cohesive  zone  model  with  rate  and  stress-state  dependencies  in 
three-dimensional  finite  element  analysis. 


1.  Introduction 

Engineering  thermoplastic  polymers  are  used  in  a  variety  of  applications.  Polycarbonate  (PC)  is  a  transparent 
engineering  polymer  noted  for  its  toughness.  Other  transparent  amorphous  polymers,  like  polymethyl  methacrylate 
(PMMA)  and  polystyrene  (PS),  exhibit  more  brittle  failure  than  PC.  Because  of  polycarbonate’s  toughness,  it  is 
useful  for  applications  such  as  protective  visors,  goggles,  helicopter  canopies,  automotive  headlight  lenses,  helmets, 
drinking  bottles,  riot  shields,  etc.  Though  PC  exhibits  greater  toughness  than  other  transparent  amorphous  polymers,  it 
is  subject  to  catastrophic  failure  in  the  presence  of  a  flaw.  In  fact,  its  impact  fracture  strength  in  the  presence  of  a  sharp 
crack  falls  into  the  range  that  is  typical  of  many  common  polymers.  Since  PC  is  used  in  protective-type  applications, 
catastrophic  failures  can  be  hazardous  to  the  individuals  involved. 

In  mode  I  loading,  polycarbonate  can  fail  in  one  of  two  ways:  ductile  or  brittle.  Hull  and  Owen  (1973)  were 
the  first  to  characterize  the  fracture  surface  morphology  of  the  two  failure  types.  Ductile  failure  can  be  described 
as  a  “tearing”  type  of  failure  with  significant  plastic  deformation  and  a  high  amount  of  energy  required  for  failure. 
When  glassy  polymers  are  loaded  under  certain  conditions  (i.e.,  plane-stress  or  in  uniaxial  compression),  they  exhibit 
shear  yielding.  Shear  yielding  is  a  bulk  process  associated  with  ductile  failure.  In  brittle  failure,  crack  growth  occurs 
rapidly  with  little  macroscopic  deformation  and  energy  dissipation.  Under  certain  tensile  conditions  (i.e.,  plane-strain 
or  faster  strain  rates),  many  polymers  fail  in  a  brittle  manner.  Brittle  failure  is  often  associated  with  crazing  in  the 
material  (Kambour  (1973),  Ishikawa  etal.  (1977),  Kramer  (1984),  and  Estevez  etal.  (2000)).  The  ductile  to  brittle 
transition  in  glassy  polymers  is  typically  viewed  as  a  competition  between  shear  yielding  and  crazing.  Crazing  is 
a  micro-mechanical  process  where  polymer  chains  align  in  the  direction  of  the  maximum  principal  stress.  Kramer 
(1984)  provides  a  comprehesive  treatise  on  the  crazing  process.  Work  performed  by  Wolstenholme  etal.  (1964), 


Ravetti  et  al.  (1975),  Broutman  and  Krishnakumar  (1976),  Yee  (1977),  and  Selden  (1987)  indicate  that  polycarbonate 
(and  some  other  glassy  polymers)  actually  have  two  distinct  energy  release  rates  associated  with  mode  I  crack  growth. 
The  lower  (or  critical)  value  is  associated  with  fast,  unstable  fracture,  while  a  higher  threshold  value  exists  that  is 
associated  with  slower,  ductile  and  potentially  stable  crack  growth.  Some  conditions  that  can  lead  to  a  transition 
from  one  failure  energy  to  the  other  include  loading  rate,  test  temperature,  plastic  zone  size,  specimen  thickness  (i.e., 
plane-stress  or  plane-strain  stress  state),  and  notch  sharpness. 

Commonly  used  approaches  for  simulating  failure  with  finite  elements  can  include  implicit  approaches  (like  the 
X-FEM  approach  introduced  by  Moes  et  al.  (1999)),  hybrid  methods  like  an  element  deletion  technique,  and  explicit 
approaches  that  use  cohesive  elements.  Implicit  approaches  utilizing  the  X-FEM  typically  enforce  a  failure  criterion 
that  results  in  element  enrichment  to  model  crack  growth.  The  method  can  be  desirable  because  exact  theoretical 
solutions  can  be  exactly  matched.  Sample  failure  criteria  could  be  critical  crack  tip  opening  displacement  (CTOD) 
or  critical  stress  intensity  factor  ( Kic )■  These  techniques  have  been  used  with  great  success,  but  can  quickly  become 
overly  complex  in  three  dimensional  analyses.  An  element  deletion  technique  was  used  by  Gearing  and  Anand  (2004) 
to  model  failure  in  glassy  polymers.  Element  deletion  is  typically  used  in  finite  strain  analyses  when  elements  can 
become  distorted  and  bulk  processes  control  failure.  Mass  and  energy  conservation  can  cause  error  accumulation 
when  using  element  deletion  techniques.  Explicit  techniques  that  utilize  cohesive  zone  elements  have  been  used  with 
considerable  success  for  modeling  failure.  The  works  of  Tijssens  et  al.  (2000),  Estevez  et  al.  (2000)  introduce  a 
rate  and  temperature  dependent  cohesive  zone  for  crazing  in  glassy  polymer  failure.  Anvari  et  al.  (2006)  introduce 
a  rate  and  stress-dependent  cohesive  zone  model  designed  for  failure  in  metals.  Models  based  on  this  work  have 
been  successfully  applied  in  2D  quasi-static  investigations  of  ductile  failure  by  Banerjee  and  Manivasagam  (2009) 
and  Scheider  et  al.  (2011). 

Caruthers  et  al.  (2004)  developed  a  thermodynamically  consistent,  nonlinear  viscoelastic  bulk  constitutive  model 
based  on  a  potential  energy  clock  (PEC)  for  modeling  polymer  deformations.  It  relies  on  the  assumption  that  all  poly¬ 
mer  relaxation  rates  depend  on  the  potential  energy  of  the  system.  Typically,  polymers  are  modeled  with  hyperelastic, 
plasticity,  or  viscoelastic  models.  The  models  can  give  similar  results  under  certain  circumstances,  but  the  physical 
differences  are  significant.  Glassy  polymers  do  exhibit  a  linear  viscoelastic  response  that  does  not  significantly  change 
with  temperature  except  that  relaxation  times  change.  Because  of  its  formulation,  the  PEC  model  predicts  mechanical 
yield  as  a  natural  consequence  of  relaxation  behavior  induced  by  loading.  Plasticity  models  predict  yield,  but  need 
a  second  mechanism  to  describe  the  glass  transition.  Hyperelastic  models  do  show  a  softening  type  of  behavior,  but 
hysteresis  effects  are  not  naturally  accounted  for.  Adolf  et  al.  (2009)  developed  a  method  of  simplifying  the  PEC 
model  to  create  a  more  usable  constitutive  model  for  general  engineering  analyses.  They  termed  the  model  the  Sim¬ 
plified  PEC  model  or  SPEC.  The  main  contributions  of  the  SPEC  model  were  to  create  a  formulation  that  was  simpler 
to  parameterize  and  implement  computationally.  The  model  was  shown  by  Adolf  et  al.  (2009)  to  successfully  capture 
the  time  and  temperature  dependent  behavior  of  polymers.  It  is  attractive  because  one  set  of  material  parameters  can 
predict  a  wide  range  of  material  behaviors  like  stress  relaxation,  yield,  volume  relaxation  after  quenching  into  glass, 
and  physical  aging  for  both  thermoplastic  and  thermoset  polymers. 

In  this  work,  a  new  cohesive  zone  model  is  developed  and  implemented  to  simulate  the  ductile  to  brittle  failure 
transition  in  polycarbonate.  The  model  is  rate  and  stress-state  dependent  and  accounts  for  both  the  high  and  low  energy 
failure  modes  associated  with  ductile  and  brittle  failure,  respectively.  The  rate  dependence  of  the  model  is  capable 
of  capturing  the  failure  mode  transition  due  to  changes  in  the  external  loading  rate.  The  stress-state  dependence  of 
the  model  is  capable  of  capturing  the  failure  mode  transition  due  to  thickness  effects,  i.e.,  plane-strain  or  plane-stress. 
For  the  first  time,  these  effects  are  coupled  in  a  way  that  a  consistent  set  of  parameters  can  be  used  to  simulate 
the  transition  from  ductile  to  brittle  failure  in  a  glassy  polymer  due  to  either  or  both  of  these  effects.  The  CZM  is 
based  on  the  simple  model  introduced  by  Tvergaard  and  Hutchinson  (1992)  but  allows  for  two  distinct  energies  per 
unit  area  of  crack  growth:  1)  a  high  energy  case  associated  with  ductile  failure  and  slow  crack  growth  and  2)  a  low 
energy  case  associated  with  brittle  failure  and  fast  crack  growth.  This  is  important  because  polycarbonate  does  not 
show  a  continuous  transition  in  failure  energy  as  loading  rate  or  stress  state  change.  Rather,  an  abrupt  transition 
occurs  from  high  to  low  energy  failure  at  a  critical  point.  The  discrete  manner  of  this  implementation  naturally  lends 
itself  to  simulation  using  explicit  dynamics.  It  is  also  important  to  recognize  that  a  robust  bulk  constitutive  model  is 
necessary  to  capture  the  failure  transition  because  in  ductile  failure,  high  amounts  of  energy  are  dissipated  through 
bulk  processes.  A  cohesive  zone  model  is  used,  in  this  analysis,  to  simply  enforce  the  appropriate  energetic  and  bulk 
response  for  failure  based  on  empirical  evidence  of  strain-rate  and  stress-state  effects. 
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First,  the  form  of  the  SPEC  model  used  for  this  study  will  be  outlined  in  Section  2.  The  new  CZM  developed  for 
this  work  will  be  presented  in  Section  3.  Section  4  will  present  the  finite  element  method  and  the  initial  boundary 
value  problem  used  in  this  work.  Section  5  will  present  the  numerical  experiments  performed  to  demonstrate  the 
problem  set-up  for  mode-I  fracture  simulations  for  the  single  edge-notch  tension  (SENT)  specimen  geometry  and  the 
results  of  those  analyses.  The  final  section  will  conclude  with  a  discussion  of  the  findings  of  this  study  and  a  summary 
of  the  contributions. 

2.  Bulk  Constitutive  Model  -  Simplified  Potential  Energy  Clock  (SPEC) 

Glassy  polymers  behave  in  a  nonlinear  viscoelastic  manner  that  is  also  temperature-dependent.  Polymer  relaxation 
rates  diminish  as  the  test  temperature  decreases  below  the  glass  transition  temperature,  Tg.  Generally,  viscoelastic  and 
viscoplastic  models  lack  the  robustness  to  describe  a  full  range  of  polymer  behaviors  such  as  stress  relaxation,  yield 
under  constant  strain  rate  loading,  volume  relaxation  after  quenching  into  glass,  and  physical  aging  (see  Caruthers 
et  al.  (2004)  and  Adolf  et  al.  (2009)). 

The  following  is  a  brief  summary  of  the  PEC  and  SPEC  models.  Readers  are  referred  to  Caruthers  et  al.  (2004) 
and  Adolf  et  al.  (2009)  for  the  full  theoretical  development  of  those  models.  The  PEC  and  SPEC  models  begin  with 
the  expansion  of  the  Helmholtz  free  energy  'P  about  an  equilibrium  state  that  would  exist  at  the  current  temperature  T 
and  density  p. 
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where  H_  is  the  Hencky  strain  measure,  which  is  a  logarithmic  strain  measure,  and  I\  is  its  first  invariant  that  is 
a  function  of  volume  only.  This  avoids  volumetric  inconsistencies  that  can  arise  with  other  strain  measures.  The 
prefactors  vPi_4  depend  on  the  current  strain  and  temperature  and  are  related  to,  though  not  necessarily  equal  to  the 
bulk  modulus  K,  the  shear  modulus  G,  the  coefficient  of  thermal  expansion  a,  and  the  specific  heat  capacity  at  a 
constant  volume  Cv,  respectively.  The  functions  /i_4  are  relaxation  functions  dependent  on  the  ‘material  time’,  t* , 
that  is  dependent  on  the  potential  energy  of  the  system,  Epot,  that  can  be  described  by  a  generalized  Williams-Landel- 
Ferry  (WLF)  equation 


(2) 


where 


log  a  =  -Ci 


rpot  _  fp°' 
ref 


C'2  +  Ep°>  -  E 


pot 

ref 


(3) 


where  Ci  is  the  first  WLF  constant  and  C,  is  related  to  the  second  WLF  constant.  The  first  Piola-Kirchhoff  stress  can 
be  calculated  by 
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where  E  is  the  Green-Lagrange  strain  measure1 .  The  Cauchy  stress  can  be  found  in  the  usual  way,  cr  =  ( p/pref)  E_-  5  ■ 

Ft,  where  F  is  the  deformation  gradient  which  can  be  broken  down  into  rotated  and  unrotated  or  rotation  and  stretch 
components  as 

F  =  RU.  (5) 

The  use  of  the  Hencky  strain  measure  requires  the  calculation  of  a  logarithmic  strain  measure  and  a  fourth-order 
tensor  for  the  Piola-Kirchhoff  stress  (see  (4)).  In  order  to  simplify  these  calculations,  the  SPEC  model  utilizes  an 
approximation  of  the  Hencky  strain  rate  by  the  unrotated  rate  of  deformation  tensor,  d 


dU 

dt 


dU 


(6) 


If  the  approximate  stress  associated  with  d  is  called  cr  ,  then  the  Cauchy  stress  can  be  calculated  by 


a_  =  —  R  ■  cr  ■  flT.  (7) 

—  Pref  —  — ^  — 

In  order  to  arrive  at  the  final  Cauchy  stress  to  be  used  in  calculations,  the  Hencky  strain  of  (1)  is  replaced  by  d. 
The  model  prefactors  vPi_4  are  assumed  to  be  constant  so  that  they  take  on  relationships  to  the  physical  quantities 
discussed  earlier.  Finally,  the  relaxation  functions  /i_4  are  replaced  with  either  f,  or  fs  which  stand  for  volumetric  and 
shear  relaxation  functions.  The  general  relaxation  functions  are  either  related  to  volumetric  quantities  (bulk  modulus, 
coefficient  of  thermal  expansion,  and  specific  volumetric  heat  capacity)  or  a  shear  quantity  (i.e.,  shear  modulus). 
Adolf  et  al.  (2009)  state  that  in  previous  parameterizations,  the  volumetric  relaxations  differed  very  little  and  could  be 
lumped  into  a  single  relaxation  function.  The  final  Cauchy  stress  is  expressed  by 
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where  (.),;  =  (.)s-(.)oo  is  the  difference  between  a  material  value  in  the  glassy  state  and  in  the  rubbery  state,  the  product 
of  the  bulk  modulus  and  the  coefficient  of  thermal  expansion  is  expressed  as  L  =  Ka,  and  e  is  the  deviatoric  part  of 

—dev 

the  integrated  rate  of  deformation  tensor.  For  completeness,  (3)  can  be  recast  as 


logfl 


CiN 

C'2+N’ 


where 


C"  =  C2  [l  +  C3aZf]  =  C2 


1  + 


TrefL', 


refaZS 


r>ref 

PrefEvd 


where  C'i  is  the  second  WLF  coefficient  and 
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is  a  tensor  quantity  that  is  distinct  from  the  scalar  Epot  used  in  (3). 
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C3  and  C4  are  clock  constants  that  describe  the  dependence  of  relaxation  times  on  volume  and  applied  deformations. 
C3  produces  an  apparent  glass  transition  temperature  with  pressure  and  C4  produces  yield.  The  volumetric  and  shear 
relaxation  functions  take  the  form  of  a  stretched  exponential  function 


fvAt)  =  exp 


(14) 


where  r  and  (3  are  the  constants  of  the  relaxation  function  that  need  to  be  determined  through  experiment. 


3.  Cohesive  Zone  Model  Development 

3.1.  Introduction  to  CZMs 

In  order  to  introduce  cohesive  zone  models,  it  is  first  beneficial  to  introduce  the  concept  of  cohesive  failure.  Zhang 
et  al.  (2007)  offer  a  clear  explanation  of  CZM  theory.  At  the  crack  tip,  materials  cannot  sustain  theoretical  infinite 
stress  values  associated  with  the  stress  singularity.  The  material  softens  in  the  crack  tip  regions  where  it  is  also  subject 
to  various  micromechanical  processes  such  as  microvoid  formation,  micro-cracking,  polymer  chain  orientation,  and/or 
crazing.  This  softening  can  be  simulated  by  a  traction-separation  law  acting  in  the  cohesive  zone  in  a  plane  along  the 
path  of  potential  crack  propagation.  The  theory  of  cohesive  zones  was  first  introduced  by  Barenblatt  (1959)  and 
Dugdale  (1960)  as  a  plastic  strip  yield  model.  Within  the  cohesive  zone,  the  material  separates  to  a  distance  Aacoh  as 
a  result  of  the  stresses  acting  near  the  crack  tip.  The  superscript  a  is  used  to  indicate  the  component  of  the  traction 
tensor  of  interest,  a  =  {«,  t\ ,  h}  for  the  normal  component  of  displacement  and  the  two  tangential  components 
and  /2,  respectively.  The  cohesive  zone  surface  sustains  a  distribution  of  tractions  o"  h,  which  are  functions  of  the 
displacement  jump  across  the  surface  Aacoh.  The  relationship  between  cHfo/]  and  Aacoh  is  defined  as  the  constitutive  law 
for  the  cohesive  zone  surface. 


3.2.  Model  Background 

The  cohesive  zone  constitutive  law  used  in  this  work  is  derived  from  the  model  presented  by  Tvergaard  and 
Hutchinson  (1992).  Figure  2  displays  the  traction-separation  behavior  of  that  model  referred  to  as  the  TH  model.  In 
the  model,  5°,  and  are  shape  parameters  that  are  bounded  by  0  <  5a  <  1  where 


Cohesive  failure  occurs 
opening,  i.e.,  A''oh  >  A"p. 


<5?  = 
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when  the  normal  cohesive  opening  is  greater  than  or  equal  to  the  peak  normal  cohesive 
The  traction  response  is  calculated  by  a  piecewise  function  as 
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Figure  1:  Generic  cohesive  zone  description  (adapted  from  Zhang  et  al.  (2007)  and  Barenblatt  (1959)).  On  the  top  left,  micro-mechanical  processes 
are  depicted  as  micro  void  formation  and  crazing.  These  are  idealized  as  a  fracture  process  zone  from  a  to  b  where  a  traction-separation  law  controls 
the  normal  opening  behavior  in  that  region. 


where  0  <  6"  <  6"  <  1  because  the  separation  value  of  the  cohesive  zone  is  normalized  by  the  maximum  separation 
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Figure  2:  Traction-separation  law  from  Tvergaard  and  Hutchinson  (1992). 


Tvergaard  and  Hutchinson  (1992)  show  that  the  work  of  separation  per  unit  area  which  is  found  by  calculating  the 
area  under  the  traction-separation  curve  is  given  by: 
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3.3.  Failure  Energy 

As  stated  earlier,  polycarbonate  (and  some  other  glassy  polymers)  actually  have  two  distinct  energy  release  rates 
associated  with  crack  growth.  The  lower  (or  critical)  value  is  associated  with  fast,  unstable  fracture,  while  a  higher 
threshold  value  exists  that  is  associated  with  slower,  ductile  and  potentially  stable  crack  growth.  In  examining  the 
transition  from  ductile  to  brittle  failure,  it  is  necessary  to  formulate  a  cohesive  zone  model  that  will  account  for  both 


6 


energy  requirements  simultaneously.  That  is,  the  cohesive  zone  model  must  require  high  amounts  of  energy  in  some 
instances,  but  also  be  capable  of  capturing  the  worst-case  (minimal  energy)  scenario  for  crack  growth. 

Figure  3  demonstrates  the  normal  direction  traction-separation  behavior  of  the  new  cohesive  zone  model  and  (18) 
displays  the  traction-separation  equation. 
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(18) 


In  the  high  energy  case,  the  normal  direction  components  are  A"  =  Amax  and  <x"  =  d-max  and  in  the  lower  energy  (or 
critical)  case.  A"  =  A m,„  and  a"  =  &mi„ .  Two  different  energy  values  are  possible  depending  on  the  fracture  mode 
(ductile  or  brittle).  In  the  tangential  directions,  frn  =  &a  =  trmin  and  A'1  =  A'2  =  Amin.  The  tangential  direction  model 
parameters  could  be  adjusted  for  a  more  general  loading  scenario,  but  the  minimum  values  of  the  normal  direction 
were  used  since  only  mode  I  loading  was  considered  in  this  work.  The  remainder  of  this  section  will  focus  on  the 
normal  opening  direction  of  the  cohesive  zone  model,  since  it  is  the  most  important  in  mode  I  loading. 


Figure  3:  Normal  traction-separation  behavior  of  new  cohesive  zone  model.  Two  distinct  traction-separation  behaviors  are  possible  depending  on 
the  expected  fracture  mode  (ductile  or  brittle). 

In  this  work,  two  factors  were  used  to  determine  which  traction-separation  curve  was  used:  normalized  cohesive 
zone  opening  rate  and  the  stress  state  of  the  surrounding  bulk  material.  The  mode  I  cohesive  zone  opening  rate  is 
directly  influenced  by  the  test  loading  rate.  The  stress  state  of  the  surrounding  material  is  directly  influenced  by  the 
specimen  thickness.  The  limits  of  the  stress  state  are  represented  as  plane-stress  and  plane-strain  scenarios;  however, 
three-dimensional  problems  will  represent  a  range  of  stress  states  somewhere  between  the  limits.  Figure  4  shows 
how  the  cohesive  energy  Ecoh  for  failure  changes  when  either  the  cohesive  opening  rate  or  the  stress  state  in  the 
surrounding  material  reach  critical  levels:  Ac  and  crRC,  respectively.  crR  is  a  stress  ratio  related  to  the  amount  of  plane 
strain  present.  It  will  be  defined  more  rigorously  in  Sections  3.5  and  3.6. 

3.4.  Rate  Dependence 

Consider  Figure  4  in  the  case  of  rate  controlled  failure.  If  the  relative  opening  rate  A  is  slow,  the  higher  energy 
traction-separation  law  will  dominate.  Alternatively,  if  the  relative  opening  rate  rate  is  fast,  the  lower  energy  traction- 
separation  law  will  dominate  the  behavior.  If  the  loading  rate  fluctuates,  the  law  is  free  to  move  between  the  two 
energy  values  and  thus,  the  two  curves  of  Figure  3. 

For  a  given  time  f,  that  is  updated  to  a  new  time  t  +  dt  over  a  time  step  of  size  dt,  the  normalized  cohesive  opening 
rate  over  dt  is  calculated  as 
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Figure  4:  Cohesive  energy  as  a  function  of  cohesive  opening  rate  or  stress  state.  Higher  opening  rates  or  higher  stress  states  (plane-strain)  result  in 
lower  energy  requirements  for  failure. 


where  A"  represents  the  cohesive  separation  at  time  t  and  A"  is  the  peak  separation  resulting  in  failure  of  the  cohesive 
zone  and  can  be  either  Am,„  or  Amax  depending  on  the  energy  required  for  failure.  A  is  divided  by  the  time  step  in 
order  to  account  for  changes  in  the  amount  of  opening  that  can  occur  because  of  changes  in  the  time  step.  That  is,  in 
a  large  time  step,  more  opening  can  occur,  while  the  rate  of  opening  is  actually  the  same.  This  is  how  the  separation 
rate  was  calculated  in  this  work  to  determine  the  controlling  traction-separation  curve.  The  model  parameter  Ac  can 
be  adjusted  to  capture  the  failure  behavior  seen  in  experiments  at  different  loading  rates. 


3.5.  Stress-State  Dependence 

The  stress-state  of  the  material  ahead  of  the  crack  front  affects  crack  propagation.  In  standard  ASTM-designated 
fracture  toughness  testing,  the  thickness  of  the  sample  must  be  large  enough  to  constrain  deformation  in  the  thickness 
direction  (predominantly  plane-strain  situation)  and  result  in  a  minimal  amount  of  energy  required  for  failure.  At  least 
one  example  of  a  stress-state  dependent  cohesive  zone  has  been  implemented  by  Anvari  et  al.  (2006),  their  work  was 
performed  on  two-dimensional  plane-strain  models.  However,  Section  5.4.1  will  show  that  three-dimensional  analysis 
is  beneficial  in  analyzing  the  plane-strain/plane-stress  condition. 

Figure  5  represents  a  planar  crack  in  an  arbitrary  domain.  Assume  the  loading  is  mode  I;  therefore,  the  loading 
direction  y'  should  be  orthogonal  to  the  crack  propagation  direction  V  and  in  the  same  direction  as  the  maximum 
principal  stress.  Therefore,  for  a  Cartesian  coordinate  system,  the  last  coordinate  direction  z'  will  be  the  direction  that 
is  orthonormal  to  the  plane  x’y’ .  The  component  of  the  stress  tensor  acting  in  the  z'  direction  is  the  component  of 
interest,  since  this  is  the  out  of  plane  stress  associated  with  the  plane  stress  or  plane  strain  condition.  It  is  desirable 
to  perform  the  cohesive  zone  stress  analysis  in  the  (V, y',z'}  space.  However,  the  stress  tensor  is  usually  expressed  in 
terms  of  the  global  coordinate  space  {x,y,  z}.  A  standard  coordinate  transformation  can  be  used  to  relate  the  global 
stress  tensor  <xy-  to  the  crack-oriented  stress  tensor  cr’km  expressed  as 


,  dx'kdx’m  n  n 

Tkm  ~  aiidxt  dxt  ~  (r:it  kjl" 


(20) 


where  /3kj  =  cos  (xJ,,  xj)  is  the  cos  of  the  x',  axis  with  respect  to  the  xj  axis. 

Since  the  in-plane  components  remain  relatively  constant,  an  appropriate  value  to  determine  the  relative  magnitude 
of  the  thickness-direction  stresses  is 
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Figure  5:  Planar  crack  in  an  arbitrary  domain. 


where  crR  is  the  bulk  stress  ratio  mentioned  in  Section  3.3.  This  ratio  is  especially  convenient  for  analysis,  since  from 
the  generalized  Hooke’s  law  in  the  plane-strain  case, 

°4  =  v(‘T'«  +  crw)  (22) 

where  v  is  Poisson’s  ratio.  Therefore, 

(Tr  =  v  (plane  -  strain)  (23) 

for  a  plane-strain  linear  elastic  problem.  By  the  definition  of  plane-stress,  cr^.  =  0;  therefore, 

< tr  =  0.  (plane  -  stress)  (24) 

It  has  been  shown  that  the  mechanical  behavior  of  polymers  depends  on  their  stress,  strain,  and  temperature 
history.  This  dependence  is  characterized  by  the  material  clock  of  the  SPEC  model.  Further,  during  explicit  dynamics 
calculations,  the  stress  (and  displacement)  field  can  become  noisy  with  spurious  high  frequency  oscillations  that  are  a 
numerical  artifact  of  crack  growth.  Therefore,  we  find  it  beneficial  to  average  crR  over  time.  A  large  number  of  time 
averaging  schemes  for  crR  are  possible.  In  this  work,  crR  was  averaged  over  a  window  of  time  as 

1  f'2 

crR  =  -  I  (Tr  dt,  (25) 

h  ~  h  Jti 

where  t\  and  h  merely  represent  a  window  of  time  in  which  to  average  crR.  The  initial  time  t\  was  taken  as  the  moment 
that  the  cohesive  traction  was  50%  of  (T„„„  from  Figure  3.  The  upper  bound  t2  was  taken  to  be  the  current  time. 


3.6.  Coupled  Stress  State  and  Rate  Effects 

The  last  step  in  formulating  the  cohesive  zone  model  is  to  couple  the  stress  state  controlled  failure  and  the  rate 
controlled  failure  into  one  parameter.  The  result  of  this  coupling  was  that  a  single  fracture  mode  would  be  predicted 
by  the  model  after  taking  into  account  the  cohesive  opening  rate  and  the  stress  state  of  the  surrounding  bulk  elements. 
Recall  that  the  critical  opening  rate  and  stress  ratio  that  result  in  the  transition  from  high  to  low  failure  energy  are  Ac 
and  crRc,  respectively.  The  fracture  mode  is  given  by  the  value  of  P,  calculated  as 


P  = 


A—  +  (1-A)^- 
Vrc  A  c 


(26) 
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where  0  <  A  <  1  is  a  separate  parameter,  determined  by  the  user,  which  allows  for  user  control  regarding  how 
much  each  critical  parameter  can  affect  the  fracture  mode.  Exclusively  rate-dependent  behavior  occurs  when  A  =  0. 
Exclusively  stress  state-dependent  behavior  occurs  when  A  =  1 .  In  order  to  investigate  the  stress-state  and  rate  coupled 
case,  A  =  0.5  in  this  work.  If  P  is  calculated  to  be  greater  than  or  equal  to  1,  it  is  assumed  that  the  low  energy  criterion 
will  control  and  fast,  brittle  failure  is  expected.  This  corresponds  to  either  a  plane  strain  stress  state  or  a  fast  loading 
rate.  Alternatively,  P  <  1  indicates  that  the  specimen  is  either  in  plane  stress  or  the  loading  rate  is  slow  and  the  high 
energy  criterion  will  control  resulting  in  slow,  ductile  failure. 

As  a  summary,  the  normal  cohesive  tractions  are  calculated  by 
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for  A''o„  <  <5, A"p 

for  <5 1  A",  <  Aljh  <  <52  A" 
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(27) 


where  A"o/;  is  the  displacement  jump  across  the  cohesive  element’s  normal  direction  (direction  of  crack  opening),  &ncgh 
is  the  actual  traction  returned  by  the  cohesive  zone,  and  &"  and  A^  are  the  peak  traction  and  cohesive  opening  values 
used  in  calculations.  Those  are  determined  by 


if  (P  >  1 ) 
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4.  Initial  Boundary  Value  Problem 

In  this  work,  the  crack  path  was  specified  a  priori.  Therefore,  the  propagation  direction  was  taken  to  be  the  same 
as  the  global  x  direction,  the  loading  direction  was  the  global  y  direction,  and  the  out  of  plane  direction  was  the  z 
direction.  In  this  orientation,  the  in-plane  stresses  are  the  crxx,  cry,,,  and  cr^  components  of  the  stress  tensor.  The 
thickness  direction  axial  stress  is  then  crzz. 

Consider  an  arbitrary  cracked  domain  SI  with  natural  and  kinematic  boundary  conditions  as  seen  in  Figure  6.  The 
boundary  of  SI  is  specified  by  T,  with  appropriate  subscripts  defining  the  surfaces  where  boundary  conditions  act  on 
SI.  The  dashed  boundary  Tc  represents  the  portion  of  the  geometry  ahead  of  the  crack  tip  where  cohesive  surfaces  are 
defined.  The  governing  equilibrium  equation  for  the  entire  body  at  an  instant  in  time  can  be  written  as 

pZ?  dSl  =  I  p'u  dSl ,  (28) 

r  n  n 

where  &  are  traction  vectors,  b  are  body  forces,  and  u  represents  the  second  time  derivative  of  the  displacement  field, 
the  acceleration.  In  general,  tractions  are  defined  as 


tr  =  it  n. 


(29) 


where  n  is  a  vector  normal  to  the  boundary  the  traction  acts  upon  which  is  chosen  pointing  outward  of  SI.  Applying 
the  definition  of  the  tractions  and  the  divergence  theorem,  the  equilibrium  equation  can  be  recast  as  a  domain  integral: 


t r+pb  —  pii  j 


dSl  =  0. 


(30) 


From  30,  the  momentum  equation  can  be  extracted: 


V_  -  cr+  pb  —  pu  =  0, 
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(31) 
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Figure  6:  Cracked  domain  with  natural  and  kinematic  boundary  conditions.  Yc  represents  the  cohesive  surface. 


where  <x  is  the  Cauchy  stress  tensor,  determined  by  (8).  The  initial  boundary  value  problem  consists  of  finding  the 
displacement  field  that  satisfies  the  momentum  equation  subject  to  initial  conditions,  natural  boundary  conditions,  and 
kinematic  boundary  conditions.  The  general  initial  conditions  are  defined  by 

u(x,  0)  =  uQ(x)  in  O 
v(x,  0)  =  y0(x)  in  Q., 

where  y  =  U  is  the  velocity.  The  natural  boundary  condition  is  given  by 

u  =  g  on  Tg,  (33) 


where  g  represents  a  known  displacement  value  on  a  surface  denoted  Tg.  The  surface  denoted  Tc  represents  the 
cohesive  surface.  The  details  of  the  traction-separation  law  of  the  cohesive  surface  model  were  presented  in  Section  3. 
On  Tc,  the  closing  traction  force  can  be  expressed  as  a  function  of  the  displacement  jump  A"  h  in  the  normal  direction 
across  the  cohesive  surface  as 

£c=^(AcJ  °n  Tc.  (34) 

Lastly,  a  kinematic  boundary  condition  is  expressed  in  terms  of  velocity  which  is  obtained  by  differentiating  the 
displacement  in  time: 

u  =  y  =  h  on  Tv  (35) 

where  h  is  a  velocity  value  or  function  of  time  that  is  known  on  the  boundary  Tv. 

Combining  all  these  equations,  the  strong  form  of  the  initial  boundary  value  problem  can  be  stated  in  the  following 
way:  given  b,  g,  h,  uQ,  and  h0,  find  u  :  flx]0,  T[— >  R  such  that 
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da_ 

pu  =  -^  +  b 
u  =  g 
v  =  h 

u(x,  0)  =  M0 
m(x,  0)  =  U  0(x) 


on  Qx]0,  T[ 
on  T?x]0,  T[ 
on  Tvx]0,  T[ 
on  Tcx]0,  T[ 
reh 
x  e  Q, 


(36) 


where  ]0,  T[  represents  the  open  interval  of  time  from  0  to  final  time  T .  The  weak  form  of  the  equation  of  motion  can 
be  obtained  by  multiplying  (31)  by  a  vector- valued  test  function  w  that  represents  the  virtual  change  in  displacement. 
Then  integrating  over  the  domain  of  the  body  and  applying  the  product  rule  and  the  divergence  theorem  to  give 
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J'  £lY  +  j'  pij_ 


w  dQ. 


11 


(37) 


j'  b  Sw  dQ.  —  J"  o~  w  dF 

Q  Tex 

f 


&  w  dr. 


rc 


where  the  last  integral  term  is  the  contribution  to  the  weak  form  from  the  cohesive  surface  arising  from  the  cohesive 
opening  Aco;,  and  rext  refers  to  the  external  boundary  of  Q,  or  all  boundaries  except  rc.  That  is, 


r  =  r,.  u  r„ 


(38) 


Introducing  the  standard  Finite  Element  discretization  gives  rise  to  the  element  matrix  problem: 

M  d  +  F‘nt  =  Ff,  (39) 

— e  —ee  e 

where  the  subscript  e  represents  the  matrix  or  vector  on  the  element  level, 

dQ.  (40) 

n, 

is  the  consistent  mass  matrix, 

F’"  =  J '  f  <r  dQ  -  J'  NrfdT 
ne  rc 

describes  the  element  internal  force  vector,  and 

Ff  =  j'  Nrbdn  +  J  if  &  dr 

ne  r„, 

describes  the  elemental  external  force  vector.  The  N_  are  matrices  containing  the  finite  element  shape  functions  and  B 
are  matrices  containing  the  spatial  derivatives  of  the  shape  functions.  The  external  tractions  arising  from  the  cohesive 
surface  on  Tc  can  be  evaluated  by  the  traction  separation  law  given  by  (27).  The  elemental  equations  are  assembled  to 
give  rise  to  the  global  matrix  problem:  given  F  :]0,  T[— >  R"e«,  find  d  :]0,  T\— >  R"'«  such  that 

Md  =  Fex,-fnt  t  e]0,  T[ 

d(0)  =  d()  (43) 

d(  0)  =  f. 

This  system  was  then  solved  using  an  explicit  time  integration  scheme  with  a  lumped  mass  matrix.  The  next  section 
will  detail  the  numerical  analyses  performed  using  the  bulk  and  cohesive  zone  constitutive  models  detailed  in  this 
section. 


(41) 


(42) 


M  =  I  pNr N 
— «  J  — '  — < « 


5.  Numerical  Analysis  and  Discussion 

5.1.  Problem  Setup 

Figure  7  shows  the  single  edge  notch  tension  (SENT)  geometry  simulated  by  finite  element  analysis.  The  SENT 
geometry  was  chosen  for  a  simple  set  up  in  which  to  perform  the  parameter  study.  In  this  work,  the  crack  path  was 
specified  a  priori  along  the  surface  rc.  Three-dimensional  simulations  were  performed  on  samples  ranging  from  1  to 
4  mm  in  thickness.  The  top  portion  of  the  specimen  represents  T,,.  A  vertical  (positive  v  direction)  velocity  ranging 
from  8000  mm/min  to  500  mm/min  was  specified  on  this  boundary.  This  loading  rate  range  gives  a  corresponding 
mean  strain  rate  range  of  13.333  s-1  to  0.833  s  .  The  bottom  surface  of  the  specimen  represents  rg.  That  surface  was 
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10  mm 


Q 

10  mm 

20  mm 

Figure  7:  SENT  domain  with  natural  and  kinematic  boundary  conditions.  The  crack  path  was  specified  a  priori  along  Tc.  A  prescribed  velocity 
was  placed  on  Tv  and  displacement  was  fixed  on  rg. 


fixed  in  the  vertical  (y )  direction.  Both  Tv  and  T^  were  fixed  in  the  x  and  z  directions  in  all  simulations.  Lastly,  the 
dotted  horizontal  line  through  the  center  of  the  specimen  denotes  the  cohesive  surface  portion  of  the  boundary  Tc. 

Using  the  SPEC  model,  bulk  deformation  depends  upon  the  thermal  and  strain  history  of  the  material.  The  thermal 
history  is  defined  from  an  initial  reference  state  at  a  reference  temperature.  In  accordance  with  Adolf  et  al.  (2009), 
the  reference  temperature  used  for  polycarbonate  is  423  K  and  is  above  the  glass  transition  temperature.  The  material 
must  be  cooled  in  the  simulation  to  the  test  temperature  (i.e.,  room  temperature).  This  occurs  (similarly  to  annealing) 
over  a  rather  large  time  frame.  As  suggested  by  Chambers  (2010),  that  time  was  taken  to  be  8500  seconds.  Therefore, 
the  first  8500  seconds  of  the  simulation  were  designated  only  for  a  thermal  cool  down  and  were  performed  in  a  quasi¬ 
static  sense  with  inertial  effects  neglected.  This  allowed  for  much  larger  time  stepping.  After  the  cool  down  was 
completed,  the  specimens  were  loaded  by  the  boundary  conditions  discussed  in  the  previous  paragraph  with  inertial 
effects  considered  in  an  explicit  dynamic  scheme.  The  SPEC  model  parameters  used  in  this  study  were  taken  from 
Adolf  et  al.  (2009)  and  are  listed  in  Table  1 . 

5.2.  High  vs.  Low  Energy  Failure  -  TH  Model  Parameter  Study 

Recall  that  the  new  CZM  implemented  in  this  work  is  derived  from  the  model  developed  by  Tvergaard  and 
Hutchinson  (1992).  If  &min  =  &max  and  Amin  =  Amax,  then  that  model  is  exactly  recovered.  It  is  then  possible  to 
perform  a  parameter  study  on  <x  and  A  to  specify  how  different  failure  energies  affect  crack  growth  in  simulations. 

As  discussed  earlier,  the  works  of  Wolstenholme  et  al.  (1964),  Ravetti  et  al.  (1975),  Broutman  and  Krishnakumar 
(1976),  Yee  (1977),  and  Selden  (1987)  provide  some  insight  into  the  amount  of  energy  required  per  unit  area  of  crack 
advance  in  PC.  All  of  the  mentioned  works,  excepting  Selden  (1987),  were  performed  on  impact  tests.  In  impact  tests, 
the  energy  per  unit  area  of  crack  growth  reported  is  in  the  range  of  1.5  x  103  to  7.0  x  103  N-m/m2  for  brittle  failures 
and  in  the  range  of  52.0  x  103  to  70.0  x  103  N-m/m2  for  ductile  failures.  Selden  (1987)  performed  failure  analysis  on 
compact  tension  specimens  and  reported  failure  energy  per  unit  area  of  crack  advance  to  be  in  the  range  of  1.5  x  103 
to  2.6  x  103  N-m/m2  for  brittle  failures  and  in  the  range  of  12.3  x  103  to  14.9  x  103  N-m/m2  for  ductile  failures.  While 
test  setup  plays  a  role  in  the  failure  energy,  the  important  conclusion  is  that  brittle  or  fast  fracture  results  in  a  minimal 
energy  per  unit  area  of  crack  growth  while  ductile  or  slow  failure  results  in  a  much  higher  energy  per  unit  area  of 
crack  growth. 

Using  these  works  as  a  baseline  for  crack  growth,  the  parameters  for  fast  fracture  were  =  50  MPa  and 
Amin  =  8  x  10-5  m.  These  values  resulted  in  a  cohesive  energy  for  failure  of  3.8  x  103  N-m/m2.  The  first  value  was 
physically  intuitive  since  50  MPa  is  near,  but  not  over,  the  yield  stress  of  polycarbonate  and  brittle  failure  is  expected 
to  occur  with  minimal  finite  strain  deformation.  The  value  of  A„„„  was  chosen  so  that  the  cohesive  energy  would  be 
in  the  expected  range  for  brittle  crack  growth.  The  subsequent  values  used  to  represent  the  high  energy  or  ductile 
failure  scenario  were  trmax  =  106  MPa  and  Amax  =  3.0  x  10-4  m.  A  parameter  study  not  shown  here  led  to  these 
values.  In  this  case,  &max  was  set  high  enough  to  force  large  strain  deformations  in  the  surrounding  bulk  material. 
It  was  expected  a  priori  that  slow  crack  growth  would  manifest  itself  by  allowing  a  finite  amount  of  crack  extension 
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Table  1 :  Parameters  values  of  the  SPEC  model 


Parameter  Value  Units 


Tref 

Koo(Tref ) 

Kg(Tref) 

Q'oa(Tref') 

ag(Tref) 

dctoo/dT 

dag/dT 

Goo{Tref) 

Gg(T„f ) 

dGoJdT 

dGgldT 

Ci 

C2 

C3 

c4 

Tv 

Pv 

Ts 

Ps 


423.15 
3.2 
4.9 
0.00058 
0.000185 
6e-7 
le-7 
1 

0.65 

0.0 

-0.2 

12 

42 

2000 

15000 

20 

0.15 

0.15 

0.38 


Kelvin 

GPa 

GPa 

ppm/°C 

ppm/°C 

ppm/°C: 

ppm/°C: 

Pa 

GPa 

GPa/°C 

MPa/°C 

°C 

°c 

°C/Pa 

,-i 


(approximately  one  to  three  elements),  then  some  time  would  pass  without  crack  growth.  Stress  would  then  build 
in  the  specimen  before  more  cohesive  failure  would  occur.  The  crack  would  continue  to  progress  in  this  start/stop 
mechanism  through  the  entire  specimen.  The  parameters  chosen  for  Amax  and  &max  represent  the  first  instances  where 
the  crack  growth  proceeded  in  an  obvious  start/stop  manner.  The  cohesive  energy  associated  with  failure  in  this  case 
with  the  parameters  mentioned  was  30.2  x  103  N-m/m2.  For  completeness,  (Si  =  0.05  and  S2  =  0.95  in  this  work. 

As  a  parameter  study,  the  values  for  &max  and  Amax  can  be  placed  into  the  TH  model  to  yield  a  slow  crack  growth 
scenario.  Alternatively,  the  parameter  values  for  &min  and  Amin  can  also  be  introduced  to  yield  a  fast  crack  growth 
scenario.  Figure  8  displays  the  instantaneous  crack  velocity  as  a  function  of  crack  position  for  the  fast  case  and  the 
slow  case.  The  total  work  per  unit  area  of  crack  growth  was  evaluated  for  the  entire  simulation  by  taking  the  area 
under  the  global  load-displacement  curve  seen  in  Figure  9.  it  was  found  to  be  3.9  x  103  N-m/m2.  Accounting  for  bulk 
elastic  deformations  increased  the  energy  required  for  crack  growth,  by  about  3%.  The  total  work  per  unit  of  crack 
growth  is  still  well  within  the  range  specified  by  the  earlier  referenced  works.  The  total  time  for  crack  growth  in  the 
slow  case  is  about  1 .0  x  10-3  s  and  about  40.0  x  10-6  s  in  the  fast  growth  case. 

Table  2  summarizes  the  findings  of  the  CZM  model  parameter  study  used  to  parameterize  the  rate  and  stress  state 
dependent  cohesive  zone  model  for  this  work. 


Table  2:  Summary  of  parameter  study  of  CZM  model  behavior  for  fast  and  slow  failure. 


Low  Energy 

High  Energy 

&  (MPa) 

50 

106 

A  (m) 

8  x  10“5 

3  x  10“4 

max  vel  (m/s) 

488 

150 

avg  vel  (m/s) 

214 

9 

Ecoh  (xlO3  N-m/m2) 

3.8 

30.2 

Etot  (xlO3  N-m/m2) 

3.9 

34.5 

Emax  (%  in  bulk) 

3 

7 

14 


Figure  8:  Fast  and  slow  crack  velocities  with  TH  model.  In  the  low  energy  case,  crack  velocity  builds  to  greater  than  1/2  the  Rayleigh  wave  speed, 
or  488  m/s.  Crack  velocity  averages  about  7.6  m/s  with  a  peak  of  22.8  m/s  in  the  middle  portion  (normalized  crack  position  between  0.1  and  0.8) 
of  the  slow  crack  growth  scenario.  Velocity  is  calculated  by  the  distance  the  crack  has  advanced  since  the  last  instance  of  crack  advance  divided  by 
the  difference  in  time  between  the  two  crack  growth  occurrences. 


Figure  9:  Load-displacement  plots  for  fast  and  slow  crack  growth  cases.  The  energy  for  failure  was  3.9  x  103  N-m/m2  in  the  fast  growth  case  and 
34.5  x  103  N-m/m2  in  the  slow  growth  case. 


5.3.  Rate-Controlled  Crack  Growth 

During  rate  dependent  simulations,  the  parameter  A  of  (26)  was  set  to  be  equal  to  zero.  Rate  effects  were  in¬ 
vestigated  using  a  thin  ( 1  mm  in  thickness)  geometry  since  thickness  effects  were  not  important.  Table  3  shows  the 
predicted  fracture  mode  with  only  rate  dependence  considered  by  the  cohesive  zone  model.  As  the  loading  rate  (units 
of  mm/min)  is  decreased,  it  is  expected  that  more  ductile  or  slow  failures  will  occur.  For  any  given  value  of  Ac, 
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this  is  the  case.  The  key  feature  of  this  model  is  the  ability  to  capture  the  energy  for  failure  required  under  different 
fracture  modes.  Conventional  wisdom  indicates  that  strain-rate  sensitive  materials  harden  with  increases  in  loading 
rate.  However,  implementation  of  this  behavior  typically  is  performed  by  constructing  a  cohesive  zone  model  that 
returns  a  higher  traction  value  at  increased  loading  rates.  This  results  in  higher  energy  values  for  failure,  which  is 
not  consistent  with  the  results  reported  earlier  by  Wolstenholme  et  al.  (1964),  Ravetti  et  al.  (1975),  Broutman  and 
Krishnakumar  (1976),  Yee  (1977),  and  Selden  (1987).  Rather,  this  work  simply  enforces  a  lower  overall  energy  for 
failure  that  corresponds  to  the  worst-case  scenario  at  high  cohesive  opening  rates  above  Ac ■  This  logic  has  the  ability 
to  capture  the  rate  dependent  fracture  mode  transition  seen  in  the  laboratory  in  PC. 


Table  3:  Rate  dependence  of  new  cohesive  zone  in  fracture  simulations. 


Loading  Rate  (mm/min) 

8000 

4000 

2000 

2.40 

X 

106 

fast 

fast 

fast 

2.45 

X 

105 

fast 

fast 

fast 

2.46 

X 

105 

fast 

fast 

slow 

2.47 

X 

105 

fast 

fast 

slow 

2.50 

X 

105 

fast 

fast 

slow 

2.51 

X 

105 

fast 

fast 

slow 

2.52 

X 

105 

fast 

slow 

slow 

2.55 

X 

105 

fast 

slow 

slow 

2.56 

X 

105 

slow 

slow 

slow 

2.58 

X 

105 

slow 

slow 

slow 

Figure  10  shows  the  crack  velocities  predicted  when  Ac  =  2.45  x  105.  For  this  value  of  Ac,  all  the  loading  rates 
resulted  in  fast  crack  growth.  Conversely,  in  Figure  11,  Ac  =  2.56  x  105  results  in  all  the  loading  rates  yielding  slow, 
stable  crack  growth.  It  appears  that  the  predicted  crack  growth  velocity  is  dependent  on  the  loading  rate  in  the  slow 
crack  growth  case.  This  behavior  is  expected  because  in  slow  crack  growth  a  few  elements  open  immediately,  then 
an  amount  of  time  passes  before  the  stress  can  build  up  enough  to  open  a  few  more  elements.  A  faster  loading  rate 
simply  results  in  less  time  to  build  the  stress  to  the  required  levels  for  failure.  No  differences  were  seen  in  the  crack 
growth  at  a  single  loading  rate,  with  different  values  of  Ac,  if  the  same  fracture  mode  was  expected.  For  example,  at  a 
loading  rate  of  4000  mm/min,  if  Ac  >  2.52  x  105,  all  crack  growths  were  slow  and  identical. 

5.4.  Stress-Controlled  Crack  Growth 
5.4.1.  Preliminaries 

As  mentioned  in  Section  1 ,  the  stress  state  of  the  material  has  a  great  influence  on  the  fracture  mode  expected. 
Adding  the  stress  state  dependence  of  (21)  is  a  major  contribution  of  this  work.  Siegmund  and  Brocks  (2000),  Tijssens 
et  al.  (2000),  Estevez  et  al.  (2000),  Gearing  and  Anand  (2004),  and  Anvari  et  al.  (2006)  are  some  authors  that  have 
used  information  from  the  stress  state  in  finite  element  simulations  of  material  failures.  With  the  exception  of  Gearing 
and  Anand  (2004)  all  of  these  works  are  limited  to  two  dimensions.  Only  Siegmund  and  Brocks  (2000)  and  Anvari 
et  al.  (2006)  use  information  from  the  bulk  field’s  stress  state  to  inform  the  cohesive  zone  model’s  behavior.  However, 
their  use  of  stress  triaxility  in  two  dimensions  misses  the  most  important  aspect  of  the  stress  state:  the  out  of  plane 
stress  (i.e.,  crzz).  Table  4  illustrates  this  point.  In  these  simulations,  two  specimens  of  different  thicknesses  (1  mm 
and  6  mm)  were  loaded  to  the  same  displacement  prior  to  cohesive  crack  growth.  The  third  column  corresponds  to 
a  1mm  thick  sample  with  a  highly  refined  mesh.  That  result  will  be  discussed  shortly.  The  stress  states  of  the  bulk 
material  were  analyzed  to  determine  the  effects  that  thickness  has  on  the  stress  state.  It  was  found  that  the  in-plane 
components  {x  and  y)  were  remarkably  consistent  even  as  thickness  changed.  However,  there  is  a  significant  change 
in  the  z  component  stress  values.  If  stress  state  is  to  be  used  in  numerical  analysis  to  determine  fracture  mode,  then 
three  dimensional  simulations  must  be  used  in  order  to  capture  the  effects  the  out-of-plane  stresses.  This  is  especially 
important  when  investigating  the  transition  from  plane  strain  to  plane  stress  or  from  fast  to  slow  failure  because  the 
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Figure  10:  Comparison  of  fast  crack  growth  velocities  at  different  loading  rates.  (Ac  =  2.45  X  105) 


Table  4:  Maximum  stress  values  in  meshes  of  two  different  thicknesses 


a-(MPa) 

6  mm 

1  mm 

lmm-fine 

O"  xx 

28 

28 

35 

CTyy 

51 

51 

51 

°~zz 

18 

5 

9 

O'  xy 

20 

20 

25 

O' xz 

10 

2.8 

3.9 

Vyz 

11.7 

4.7 

7 

stress  state  in  a  three  dimensional  specimen  is  a  combination  of  plane  strain  in  the  interior  and  plane  stress  at  the 
exterior  of  the  geometry. 

Figures  12  and  13  show  the  effect  of  thickness  on  crR  in  simulations.  The  ratio  crR  calculated  from  (21)  is  a  natural 
way  to  determine  how  “close”  a  portion  of  the  geometry  is  to  a  plane  stress  or  a  plane  strain  stress  state.  The  figure 
displays  the  plane  in  which  a  crack  would  propagate.  The  black  line  represents  the  crack  front  with  the  arrows  pointing 
in  the  direction  of  crack  growth.  The  area  immediately  ahead  of  the  crack  tip  experiences  a  maximum  crR  value  of 
0.22  in  the  6  mm  sample;  the  1  mm  sample  exhibits  a  maximum  crR  of  0.062.  This  shows  that  there  is  a  quantitative 
difference  between  crR  in  specimens  of  different  sample  thicknesses.  These  meshes  are  identical  with  the  exception 
that  Figure  12  is  six  times  the  thickness  of  Figure  13.  The  element  sizes  are  all  identical.  This  results  in  a  rather  coarse 
mesh  for  the  1  mm  geometry. 

To  investigate  the  effect  that  mesh  refinement  has  on  crR,  the  1  mm  thick  sample  was  refined  to  an  element  height 
three  times  smaller  than  that  seen  here.  Therefore,  there  were  18  elements  through  the  thickness  instead  of  6  as  shown 
in  the  figure.  Figure  14  shows  the  crR  response  of  the  finer  1  mm  mesh.  The  crR  scale  is  the  same  as  the  6  mm  mesh 
from  earlier.  Notice  that  there  is  an  obvious  quantitative  difference  between  the  crR  values  in  each  geometry.  The 
peak  crR  in  the  fine  mesh  is  about  0.12.  These  values  correspond  to  a  percent  difference  of  about  83%.  The  mesh 
refinement  leads  to  a  slight  change  in  the  stress  tensor  as  seen  in  Table  4.  Additionally,  consider  the  stress  criterion 
of  a  critical  hydrostatic  stress  being  used  to  determine  if  plane  stress  or  plane  strain  controls.  Recall  that  hydrostatic 
stress  is  defined  as 

CrH  =  ^  (O'™  +  0-yy  +  Czz)-  (44) 
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Figure  1 1:  Comparison  of  slow  crack  growth  velocities  at  different  loading  rates.  (Ac  =  2.56  X  105) 

From  Table  4,  the  crH  value  for  each  mesh  is  97,  84,  and  95  MPa  respectively.  There  appears  to  be  a  difference 
between  the  1mm  and  6mm  thick  geometries;  that  difference  is  about  15%  initially  for  the  same  level  of  discretization. 
However,  when  the  1mm  mesh  is  refined,  that  difference  quickly  disappears. 
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Figure  12:  Stress  ratio  ( ctr )  in  6  mm  sample. 


In  the  previous  example,  crR  was  evaluated  by  analyzing  the  bulk  elements  that  border  cohesive  elements.  It  is 
necessary  to  inform  the  cohesive  zone  elements  with  crR  from  the  surrounding  bulk  elements  to  perform  this  type 
of  analysis.  As  an  example,  consider  a  cohesive  element  that  is  bordered  by  two  adjacent  bulk  elements  in  a  three- 
dimensional  mesh  as  seen  in  Figure  15.  During  the  cohesive  surface  initialization,  the  bulk  element  global  identifi¬ 
cation  numbers  were  determined  by  searching  the  mesh  for  shared  nodes.  For  the  given  example,  this  process  would 
determine  that  cohesive  element  #3  is  attached  to  bulk  elements  #1  and  #2.  When  cohesive  element  #3  is  instructed  to 
provide  its  traction  value  for  the  internal  force  vector,  the  current  stress  of  state  of  elements  #2  and  #1  are  processed 
to  return  a  creR  value  for  each  element.  Then,  both  creR  values  are  averaged  to  return  the  average  bulk  stress  ratio  crR 
actually  experienced  by  the  cohesive  surface  element.  This  value  is  then  passed  to  the  cohesive  surface  model  routine 
to  determine  whether  a  high  or  low  energy  threshold  should  be  used  in  the  calculation. 

The  method  implemented  to  return  crR  is  only  valid  if  no  topological  mesh  changes  occur  during  the  simulation. 
That  is  the  case  for  the  simulations  run  in  this  work.  However,  if  cohesive  surface  elements  were  to  be  dynamically 
inserted  into  the  mesh,  the  search  for  bordering  elements  would  have  to  be  performed  just  before  calculating  the 
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Figure  13:  Stress  ratio  (ctr)  for  coarse  mesh  in  1mm  sample. 
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Figure  14:  Stress  ratio  ( ctr )  for  fine  mesh  in  1mm  sample. 


Figure  15:  Cohesive  surface  element  inserted  between  two  bulk  finite  elements.  The  cohesive  surface  relies  on  stress  state  information  from  the 
surrounding  bulk  elements  to  determine  the  failure  energy  required. 


internal  force  vector,  not  during  the  initialize  step.  This  would  result  in  significant  reductions  in  a  finite  element 
code’s  performance. 
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5.4.2.  Stress  Dependent  Cohesive  Zone  Results 

In  order  to  elucidate  the  stress  state  dependence  of  the  model,  the  crR  term  was  isolated  by  placing  A  =  1  in  (26). 
Simulations  were  performed  in  the  same  manner  as  in  the  rate  dependent  case.  In  this  parameter  study,  o-RC  was  varied 
as  was  the  specimen  thickness.  This  resulted  in  determining  of  crRc  values  where  each  specimen  thickness  transitions 
from  slow  to  fast  fracture. 


Table  5:  Stress  dependence  of  new  cohesive  zone  in  fracture  simulations 


O' RC 

1  mm 

2  mm 

3  mm 

0.007 

fast 

fast 

fast 

0.008 

trans 

fast 

fast 

0.009 

trans 

trans 

fast 

0.010 

slow 

trans 

fast 

0.011 

slow 

trans 

fast 

0.012 

slow 

trans 

trans 

0.013 

slow 

slow 

trans 

0.014 

slow 

slow 

slow 

0.015 

slow 

slow 

slow 

0.016 

slow 

slow 

slow 

Table  5  shows  the  results  of  the  crRC  and  thickness  parameter  study.  The  fracture  modes  “fast”  and  “slow”  are 
defined  in  the  same  way  as  the  previous  section.  However,  now  there  is  a  separate  fracture  mode  termed  “trans”  for 
transition.  The  transition  fracture  mode  indicates  that  one  particular  type  of  fracture  mode  did  not  dominate  the  entire 
fracture.  For  example,  in  some  cases  a  crack  would  progress  slowly  for  some  time,  then  quickly  through  the  remainder 
of  the  specimen.  If  the  slow  crack  growth  portion  was  less  than  1/2  of  the  specimen  fracture  area,  fast  crack  growth 
was  said  to  be  the  controlling  mode.  In  simulations,  if  fast  fracture  would  occur  after  1/2  of  the  total  area  of  crack 
growth  had  occurred,  then  the  failure  was  called  a  transition  failure.  The  one  exception  to  this  occurred  when  the 
last  few  elements  near  the  end  wall  would  quickly  fail.  This  was  called  slow  fracture  because  the  speed  of  the  crack 
advance  was  overwhelmingly  slow  and  the  failure  energy  was  obviously  high  energy,  but  only  the  last  3-6  rows  of 
elements  would  quickly  open.  This  slow  behavior  was  explained  earlier  and  can  clearly  be  seen  in  Figure  8 

Figure  16  shows  the  simulated  fracture  velocities  of  a  fast,  transition,  and  slow  failure  in  a  2  mm  thick  sample 
using  the  stress  controlled  cohesive  zone  model.  The  stress  controlled  fast  and  transition  failures  show  a  period  of 
slower  velocity  crack  growth  at  the  beginning  of  the  fracture.  This  period  seems  to  correspond  very  well  to  zone  I 
of  PC’s  fracture  surface  morphology  characterized  by  Hull  and  Owen  (1973).  During  that  time,  stress  builds  and  the 
crack  progresses  slowly  in  a  relative  thumbnail  shape.  Eventually,  in  the  transition  or  fast  cases,  the  crack  will  begin 
to  progress  quickly.  In  the  slow  cases,  it  continues  to  slowly  progress. 

An  example  of  the  thumbnail  shaped  crack  can  be  seen  in  Figure  17,  which  shows  a  2  mm  thick  specimen  near  the 
end  of  zone  I,  just  before  the  crack  progresses  quickly.  The  top  half  of  the  specimen  is  removed  so  that  the  cohesive 
zone  elements  along  the  plane  of  crack  growth  may  be  viewed.  The  bottom  half  of  the  specimen  is  colored  by  the 
maximum  principal  stress  values  per  element.  The  cohesive  elements  are  colored  gray  and  are  set  to  disappear  after 
cohesive  death  has  occurred.  This  enables  the  visualization  of  the  crack  tip  as  it  advances. 

If  a  single  crRC  value  is  considered,  like  0.011,  it  is  easy  to  see  that  a  transition  in  fracture  mode  occurs  as  the 
specimen  thickness  is  increased.  In  the  next  section,  stress  state  and  rate  effects  will  be  coupled  to  simulate  a  fracture 
mode  transition  with  both  effects  considered  simultaneously. 

5.5.  Coupled  Rate  and  Stress  Effects 

The  coupling  of  rate  and  stress  state  effects  by  (26)  is  the  key  feature  of  this  model.  The  final  parameters  investi¬ 
gated  in  this  study  were  A  =  0.5,  crRC  =  0.0065,  and  Ac  =  2.50  x  105.  Table  6  shows  the  results  for  fracture  mode  for 
the  material  parameters  given. 

The  parameters  crRC  and  Ac  were  chosen  with  time  constraints  in  mind.  That  is,  crRC  was  chosen  lower  than  the 
expected  value,  based  on  the  study  from  Section  5.4.2.  This  allowed  for  the  transition  from  ductile  to  brittle  failure 
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Figure  16:  Simulated  fracture  velocities  with  stress  controlled  cohesive  zone.  The  crack  velocity  transitions  from  fast  to  slow  as  ctr  is  adjusted  to 
require  a  stress  state  approaching  plane  strain  for  low  energy  failure. 


Figure  17:  Zone  I  and  a  thumbnail-shaped  crack  front  in  a  finite  element  simulation. 


to  occur  at  a  lower  value  of  thickness.  Therefore,  fewer  total  elements  were  needed  in  simulations.  For  example,  all 
other  things  remaining  equal,  a  mesh  that  is  1  mm  in  thickness  will  have  three  times  fewer  elements  than  a  mesh  that 
is  3  mm  in  thickness.  This  explains  the  very  minor  changes  in  the  fracture  mode  results  at  thicknesses  greater  than 
2  mm.  Likewise,  Ac  was  chosen  in  such  a  way  as  to  minimize  simulation  time.  This  resulted  in  an  arbitrarily  large 
value  for  the  parameter,  but  allowed  for  a  fracture  mode  transition  to  occur  at  higher  overall  loading  rates.  In  explicit 
dynamics,  the  time  step  size  is  limited  for  solution  stability.  Therefore  a  simulation  run  with  an  8000  mm/min  loading 
rate  actually  takes  about  eight  times  shorter  to  run  than  a  simulation  with  a  1000  mm/min  loading  rate.  To  further 
appreciate  the  time  constraints  mentioned  above,  a  1  mm  thick  sample  loaded  at  8000  mm/min  had  a  run  time  of  just 
over  7  hours  when  run  on  16  processors  of  the  MPI  cluster  used  in  this  work. 

In  summary,  the  data  of  Table  6  display  unprecedented  results.  Consider  the  fast  fracture  mode  of  the  1  mm  thick 
geometry  loaded  at  8000  mm/min.  Due  to  thickness,  it  would  be  expected  that  a  1  mm  sample  is  very  nearly  in  a 
plane-stress  stress-state  and  slow  fracture  would  be  expected.  However,  due  to  the  fast  loading  rate  dominating  in  the 
controlling  equation  (26),  the  sample  fails  in  a  brittle  manner.  This  is  consistent  with  the  behavior  seen  in  thin  samples 
of  PC.  As  the  loading  rate  is  decreased  (e.g.,  6000  mm/min),  the  fracture  mode  transitions  into  the  slow  regime.  This 
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Table  6:  Fracture  mode  transition  for  coupled  problem. 


Loading  Rate 

1  mm 

2  mm 

3  mm 

4  mm 

8000 

fast 

fast 

fast 

fast 

7800 

fast 

fast 

fast 

fast 

7600 

fast 

fast 

fast 

fast 

7500 

slow 

fast 

fast 

fast 

7000 

slow 

fast 

fast 

fast 

6000 

slow 

fast 

fast 

fast 

5000 

slow 

fast 

fast 

trans 

4000 

slow 

trans 

trans 

trans 

2000 

slow 

trans 

trans 

trans 

1000 

slow 

trans 

slow 

trans 

500 

slow 

slow 

slow 

slow 

allows  for  the  stress-state  in  the  material  ahead  of  the  crack  tip  to  control  the  failure  mode.  Then  as  the  thickness  is 
increased,  but  loading  rate  is  held  constant,  the  fracture  mode  transitions  back  to  fast  fracture.  This  occurs  because  the 
stress-state  has  been  increased  toward  a  plane-strain  stress-state.  This  is  the  first  instance  of  using  numerical  analysis 
to  simulate  a  fracture  mode  transition  due  to  rate  and  thickness  effects  with  a  single  set  of  material  input  parameters. 


6.  Conclusions 

In  this  work,  the  SPEC  model  has  been  applied  to  mode  I  fracture  problems  in  thermoplastic  glassy  polymers. 
The  model  is  a  continuum  model;  therefore,  it  is  not  designed  to  capture  the  different  micromechanical  processes, 
such  as  crazing,  that  influence  the  deformation  behavior  of  the  material.  However,  the  new  cohesive  zone  model 
was  implemented  to  phenomenologically  capture  the  micromechanical  processes  in  the  fracture  process  zone  that  the 
SPEC  model  was  not  designed  to  capture. 

A  single  set  of  input  parameters  have  been  used  to  simulate  the  transition  from  slow  to  fast  fracture  as  both 
thickness  and  loading  rate  change.  Three-dimensional  explicit  dynamic  finite  element  analysis  was  performed  and 
a  new  cohesive  zone  model  was  implemented  that  is  rate  and  stress  state  dependent.  Therefore,  loading  rate  and 
thickness  inform  the  traction-separation  behavior.  The  cohesive  zone  model’s  behavior  accounts  for  both  the  low 
energy  critical  scenario  of  fast,  brittle  crack  growth  and  the  high  energy  scenario  of  slow,  ductile  crack  growth.  While 
there  are  some  rate  dependent  and/or  stress  dependent  cohesive  models,  these  models  either  are  not  applicable  to 
polymer  behavior,  or  are  not  capable  of  capturing  a  fracture  mode  transition. 

In  this  work,  the  cohesive  traction  remains  constant  in  areas  of  either  ductile  or  brittle  failure  and  reduces  from 
high  to  low  as  loading  rate  is  increased.  This  results  in  less  bulk  yielding  at  high  loading  rates  and  thus  more  fast, 
brittle  failures.  The  rate  dependence  of  the  cohesive  model  is  incorporated  by  determining  the  cohesive  opening  rate. 
The  model  uses  a  critical  cohesive  opening  rate  Ac  as  the  threshold  opening  rate.  Below  Ac,  the  high  energy  regime 
controls,  resulting  in  slow,  ductile  failure.  Above  Ac,  the  low  energy  regime  controls,  resulting  in  fast,  brittle  failure. 

Three-dimensional  analysis  is  of  primary  importance  when  performing  stress  analysis  aimed  at  the  ductile  to  brittle 
transition.  Plane-strain  results  in  fast,  brittle  failure  and  plane-stress  results  in  slow,  ductile  failure.  Two-dimensional 
analysis  idealizes  stress  in  the  third  dimension.  Previous  stress  state  dependent  cohesive  zone  models  are  limited  by 
two-dimensional  analysis. 

Here  we  show  that  regardless  of  thickness,  in-plane  stresses  do  not  vary  significantly.  However,  the  stress  in  the 
thickness  direction  varies  significantly.  A  stress  ratio  crR  was  proposed  that  is  able  to  quantify  the  amount  of  plane- 
strain  or  plane-stress  that  exists  at  a  point.  A  critical  crRC  was  utilized  above  which  the  low  energy  regime  controlled 
crack  growth  and  above  which  the  high  energy  regime  controlled  crack  growth. 

The  stress  and  rate  dependencies  were  coupled  in  a  simple,  linear  way.  Without  changing  input  parameters,  the 
model  shows  a  fast,  brittle  failure  in  a  thin  sample  thickness  when  loaded  at  a  high  loading  rate.  As  the  loading  rate 
is  reduced,  the  model  transitions  to  slow,  ductile  crack  growth.  As  the  sample  thickness  is  increased,  while  holding 
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loading  rate  constant  at  the  lower  value,  the  model  transitions  back  to  the  fast  crack  growth  regime.  This  behavior 
is  consistent  with  the  behavior  seen  in  the  experimental  analysis.  The  model  can  easily  be  used  and  compared  with 
experimental  results  to  construct  the  full  fracture  mode  transition  behavior  based  on  loading  rate  and  thickness  effects. 
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